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Abstract 



Nuclear reactions govern major aspects of the chemical evolution of galaxies and stars. Analytic study 
of the reaction rates and reaction probability integrals is attempted here. Exact expressions for the reaction 
rates and reaction probability integrals for nuclear reactions in the cases of nonresonant, modified nonres- 
onant, screened nonresonant and resonant cases are given. These are expressed in terms of H-functions, 
G-functions and in computable series forms. Computational aspects arc also discussed. 
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I. Introduction 



Nuclear reactions govern major aspects of the chemical evolution of the universe or, at least, its building 
blocks: galaxies and stars. A proper understanding of the nuclear reactions that are going on in hot 
cosmic plasmas, and those in the laboratories as well, requires a sound theory of nuclear-reaction dynamics 
(Brown and Jarmie 1 ). The nuclear reaction rate is the central quantity connecting the theoretical models of 
universal, galactic and stellar evolution and the nucleosynthesizing nuclear reactions, which can be studied 
to some extent in nuclear-physics laboratories. Compilations of reaction rates and uncertainties (analytic 
expressions and tabulated temperature-dependent values) and astrophysical 5-factors (analytic expressions 
and tabulated energy-dependent values) for charged-particle induced reactions are available on the World 



Wide Web (http://csa5.lbl.gov/chu/astro.html). The rate of reacting particles i and j, in the case of 



nonrclativistic nuclear reactions taking place in a nondegenerate environment, is usually expressed as 



7r/i J \kT 

= TliTlj < av > (1-1) 

where rii and rij denote the particle number densities of the reacting particles i and j, fi = . is the 

reduced mass of the reacting particles, T is the temperature, k is the Boltzmann constant, a(E) is the 
reaction cross section and v = {2E/p)^ is the relative velocity. Thus < av > is the reaction probability 
integral, that is, the probability per unit time that two particles, confined to a unit volume, will react with 
each other. The basic assumptions made in (1.1) are that the reacting particles i and j have isotropic 
Maxwell-Boltzmann kinetic-energy distributions and that the kinematic of the reaction can be treated in the 
center-of-mass system, see Fowler 2 , Mathai and Haubold 3 , and Imshennik 4 . 

For nonresonant nuclear reactions between nuclei of charges Zi and Zj at low energies (below the Coulomb 
barrier), the reaction cross section has the form 

a ( E ) = ^El c -^v(E) 
v ' E 



ith 



2/ UE 1 / 2 



where rj(E) is the Sommerfeld parameter, fi is Planck's quantum of action and e is the quantum of electric 
charge. For a slowly varying cross-section factor S(E) we may expand 

S(E) = S(0 ) + ^« £ + I^^ 



dE 2 dE 2 

^(fcr)-^ v\ 

x / y v e- y &- zy ^dy 



Thus 

< av > 



where 

E ( n \ I z l z J e 2 

Then the integral to be evaluated is 

r°° _i 

N v (z) = / y v e- y e- z y 2 dy. (1.2) 



o 



We will consider a general integral of the following form: 
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Nonresonant case: 

h = / y u e- ay e- zy ~"dy, a > 0, z > 0, p > 0. (1.3) 
Jo 

It may be a stringent assumption to consider a thermonuclear fusion plasma as being in thermodynamic 
equilibrium. If there is a cut-off of the high energy tail of the Maxwcll-Boltzmann distribution then (1.2) 
becomes the following: 

Nonresonant case with high energy cut-off: 

r d 

T 2 = / y v e- ay e~ zy , d<oo, a > 0, z > 0, p > 0. (1.4) 
Jo 

We consider also ad hoc modifications of the Maxwell-Boltzmann distribution for the evaluation of the 
nonresonant thermonuclear reaction rate, which acts as a depletion of the tail of the distribution function 
(see Kaniadakis et al. 5 In this case the general integral to be evaluated is the following: 

Nonresonant case with depleted tail: 

/>oo 

h= y"e- ay e- byS e- zy ~''dy, 5 > 0, z > 0, p > 0, a > 0. b > 0. (1.5) 
Jo 

Taking into account the electron screening effects for the reacting particles, Haubold and Mathai 6 
consider the case of the screened nonresonant nuclear reaction rate (see Lapenta and Quarati 7 ). In this case 



x / y v e- y e-< y+ —) 2 dy. (1.6) 
Jo 

In this case we consider the following general integral: 
Screened nonresonant case: 

/>OC 

h = / y v e- ay e- z{y+t) ~ P dy, t > 0, p > 0, z > 0, a > 0. (1.7) 
Jo 

When the nuclear cross section <r(E) in (1.1) has a broad single resonance it can be expressed via the 
parametrized Breit-Wigner formula and then we have, see also Haubold and Mathai 8 , 



< av > = (2tt) ^ 



£ ZiZ 3 e 2 RowT k iD 



i + arM 2 / 



00 p-o-y-qy' 



(i ri ) 2 y (b-y) 2 + g [ 



■dy (1.8) 



where 



b — E r — -r r!, 



fcT(l + (I ri ) 2 ) 4 
g=\(T + E r T l ), q = q^l+(^T 1 



Thus the general integral to be evaluated in this case is of the following form: 
Resonant case: 

h = / 77 w— -dy- (1.9) 

Jo (b - V? + 9 

In the resonant case also we can consider a modification of the Maxwell-Boltzmann distribution which 
results in a depletion of the tail. Then the general integral to be evaluated is the following: 

Resonant case with the depleted tail: 

J 6 = / — ( xT—^Ay. (1.10) 

Jo (c - vr + g 2 



The aim of the present article is to give new exact analytic representations of the integrals in (1.2) to 
(1.10). Additional representations are available from Mathai and Haubold 3 . 

II. Reduction formulae for the reaction probability integrals 

Writing 



we have 
Now consider 



r(oo) _ j 

1 2 — 1\. 



< e -ay e -by s e -zy » dy 



/•OC 

= / 2A 
Jo 

= ^ — — / y v + Sm e - av e- zv "dy 

m =0 m ' ^° 

°° ( — h\ m , n 

= T,- J T I t\v + 5rn,a,z,p) (2.1) 



m! 

m—0 



where (a)„ denotes the Pochhammcr symbol, 

(a) n = a(a + l)...(a + n - 1), (a) = 1, a ^ 0. 

Note that 



y u e -ay e -z(y+t) » dy 

/OO 
(u - ty 

±fn ™ ! Ju=t 



t 

at 

m=0 

cat E tm (- m ' a < z < p) - T 2 ] (-™< z ' ")] • ( 2 - 2 ) 

m=0 



For simplifying 25 and Iq we will use the identity 



Jo 



(c - y) 2 + g 2 
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and rewrite the single integral as a double integral. That is, 



Now we expand 



roc poo 

Jx=o Jy=0 



ml 

m=0 

°° 2m /9m\ r 2ni- mi 

E E CY-^ m+mi —^- m v mi 



ra— mi— 



mil m\ 



The integral over x gives 



Substituting back we have 



/•OO 

/ a^c^dx = (. 9 2 )-(" l+1 )m!. (2.3) 

>/x=0 



_ 1 y, ^ (2m)!c 2m - mi (-l) m+mi 
5 ~ ^2 2^ 2^ mi !(2m-mi)! (g 2 ) m 

y m=0mi=0 v ; va ; 

/■oo 

x / y I/+mi e- ay e-^"' > dy 
./o 

oo 2m 



y m=0mi=0 v 7 \ H / 



Thus it is seen that all the integrals 1\ to I e can be reduced to the integral I^"\v, a, z, p) for two different 
situations of non-negative as well as negative v. We will evaluate in the next section by using a statistical 
technique. 

III. Evaluation of the integral 

In general, integrals I\ to Iq are quite difficult to evaluate analytically. Here we will use a statistical 
technique. We will evaluate the density of a product of two independently distributed real scalar random 
variables by using two different methods, one procedure leading to the integral that we want to evaluate 
and the other procedure leading to a representation in terms of a known function. Then appealing to the 
uniqueness of the density we have the integral evaluated in terms of a known special function. Let x and y 
be real scalar random variables having the densities 



and 



ae- ax , < x < d 
0, elsewhere 



My) : 

where c\ and c-i are normalizing factors such that 



c 2 y v e~ zy , < y < oo 
0, elsewhere 



r-d r-oo 

/ h(x)dx = 1 and / f 2 (y)dy = 1. 
Jx=0 Jy=0 
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Since the variables arc assumed to be independently distributed, the joint density of x and y is the product 
of f\{x) and f2(y)- Let us transform x and y to u = xy and v = x. Then the joint density of u and v, 
denoted by g(u, v), and the marginal density of u, denoted by gi(u), are given by 



v n .—v—\ n —av ~ — cv 



e e 



c = ZU 1 ' 



and 



Hence we have 



9i{u) = c lC2 u v / v- , - 1 e- av e- co "dv, c = zu». 



/ v- v ~h 
Jo 



cic 2 



gi(u), c = zu p . 



(3.1) 



Let us look at the (s — l)-th moment of u. Due to independence 

^(u- 1 ) = [Eix'-^W- 1 )]. 



But 



r d 

Eix 3 - 1 ) =ci / x'^e-^dx 
Jo 

(-l) m (ad) m d s 



m=0 



s + m 



for cZ < 00 



= cia- s r(s),9?(s) > 0, for d = 00 
where SFf(-) denotes the real par of (•), and 



/>oo 

Eiy 3 - 1 ) = c 2 / j,"+- 1 e-^di/ 
Jo 

, for H x ;/' + s) > 0. 



C2 -i±£-p y + s 

— z p r 



Taking the inverse Mellin transform of E(u s 1 ) we have 



, > cic 2 _, 
gi(u) = z p 



2m J L 



E 

m=0 
d s 



{-l) m {ad) r 



L s + m 



Z~pT 



V + s 



u s ds 



(3.2) 



where L is a suitable contour and i = \/—T. This contour integral can be written as an H-function, see for 
example Mathai and Saxena 9 . That is, 



1 



s + m 



-Z~~pT 



p 



u- s ds = Hl§ 



1 

UZP 


(m+1,1) 


d 


(*m>.(^). 



Substituting (3.3) in (3.2) and then comparing with (3.1) we have 



v - v -l e -av e -zu"v " 



dv = {—v — 1, a, zu p , p) 



(3.3) 



(3.4) 



z PU 



E 



{-ad) 1 



p ' m\ 

F m=0 



x H 



1 

UZP 


(m+1,1) 


d 


("MM). 



1.2 



z pu - 20 

n 0.2 



, for d < 00 



1 




uazp 


(°.i>.(*.*). 



for d = 00 



(3.5) 
(3.6) 



where $l(v) > 0,SR(o) > 0, 5R(z) > 0,3fc(p) > 0,0 < u < oo. Note that (3.4) , (3.5) and (3.6) give an 
lf\v,a, z, p) with v negative. If the integral for a positive v is required then we proceed as follows: Take 
fi(x) = c 3 x u e~ ax and f2(y) — c±c~ zx " , where C3 and C4 are the new normalizing constants, and then proceed 
as before. Then we end up with the following representations. 



f 

Jo 



1 e - av e-* v "dv = l!p(v-l,a,z,p) 



p , 



{-ad) r 
m! 



x if J 2 



2.0 
2 



~d 



(1/ + m + 1, 1) 
(z/ + m,l), 



- a U 2 ' 

- —H , 2 



azp 



(".i),(o,i) 



(<••*) 

for = 00 



for cZ < 00 



(3.7) 



(3.8) 
(3.9) 



where 5R(z/) > 0, SR(a) > 0, 5R(,z) > 0, 5R(p) > 0. Thus (3.7), (3.8) and (3.9) cover the case of a positive v 



in 7 2 (•)• When p is real and rational then the H-functions appearing in (3.3) to (3.9) can be reduced to 
G-functions by using the multiplication formula for gamma functions, namely, 



T{mz) = (27r)^^m r ' 









m — 1 \ 








...r (zh 




, m 




m / 




m ) 





(3.10) 



For the theory and applications of G-functions see for example Mathai 1 0. For p = 21, m,n — 1,2,... 
reduction to the G- function is available from Mathai and Haubold 3 . 

The parameters of interest in nuclear astrophysics are I^"\v,a,z,p) for a = 1, z > 0, p = \. In this 
case computable representations will be discussed in the next section. 



IV. Computable series representations of the reaction rate integrals 

Let us start with (3.9) for p — \. Then the H-function to be evaluated is the following: 

H ^ { ' ] = ^nj r(^ + S )r(2 S )(az 2 )-M S 

/ 2 \ — s 

ds 



1 1 

V^2^ri 



r( s )r( s + -|r(, + s ) 



az 



(4.1) 



by expanding T(2s) with the help of (3.10). Note that (4.1) is a G-function of the type Gq'®(-), see for 

A 

2 1 



example Mathai 10 . Observe that for v ^ ± ^, A = 0, 1, ... all the poles of the integrand are simple. Then 



(4.1) generates a series corresponding to each gamma in the integrand. Corresponding to T(s) the poles are 
at s = — n, n = 0, 1, ... and the corresponding residue is 



lim 



(-iy 



-,)■ -'<) 



0? 



But 



and 



T{y-n) = 



(-i)-rQ) 

(-1) T T» 

(l-f)n ' 
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Hence the sum of the residues gives 



(-1)" 2f- 



n=0 \2)n y 1 U ' r 



T(u) F 2 



1 az z 



(4.2) 



where 0-^2 (■) is a hypergeometric series which is convergent for all a and b. Thus (3.9) is a linear function 
of 3 such series of o-PYs for v ^ n = 0, 1, .... If ^ is an integer or half-integer then we have one set of 
poles of order 2 each and one set of order one each. Techniques for handling the integral when the integrand 
has higher order poles are given in Mathai 10 . Series representations of (4.1) for all cases of v are given in 
Mathai and Haubold 3 . 



H = H*$ 



Now let us examine (3.8) for p = i. In this case the H-function to be evaluated is given by 

(i/+m,l),(0,2) 



ds. 



2m J v + m + s v ' \ d J 

-I — 

2tri J v + rr 



m + s 



T(s)T s 



(4.3) 



At s = — v — m there is a pole of order one if v ^ A = 0, 1, otherwise it will be a pole of order 2 at 
this point. When it is a pole of order one the residue is given by 



2^ 



T(-u-m)T y-v-m 
1 T{-v)T(-v- 



2 \ m-\-v 



d 

2 \ m-\-v 



~ 2^(v + l) m ( v+ \) m \ d 
Hence corresponding to this residue the term in (3.8) is the following: 



V^ (v + l)m («/+£), 



-{-az z ) 



2\m 



^r ( -„r(-, + i 



0^2 



; v + 1, v + -; -az 2 



which is convergent. 

At s = —n, n = 0, 1, ... the integrand in (4.3) has poles of order one when v ±±\, A = 0,1,... The 



sum of the residues here is given by 

77^ S 



(-1)" 



2v7r ^-^ + m — n n\ 

v n=0 



r -n 



1 00 

= 5E 



1 



1 



2\ ™ 



2^ j/ + m-n(i) n! 

n=0 V2/ n 

The term corresponding to this in (3.8) is the following double series: 

{-ad) m 1 1 (z 2 /d) r 



(4.4) 



OO OO 



m— n— 



m 



! „ + TO - n (§) n 



n 



8 



By Horn's theorem on convergence, see for example Srivastava and Karlsson 11 this double series is convergent 
for all a, z and d. The third term of (3.8) as well as all terms for other cases of v can be seen to give convergent 
series for all < d < oo, a > 0, z > 0. Similar arguments hold good for (3.5) and (3.6) also. 

Let us examine I3. We have expressed ^3 in terms of in (2.1). That is, 



m=0 



For p = i a typical term in this behaves like a 0-F2 given in (4.2). A typical term will be a constant 
multiple of a double series of the form 



00 00 



EE 

711—0 n—0 



(-by 



(-*)' 



(l)„( 1 - I/ - (5m )™ n! 



(4.4) 



This by Horn's theorem is convergent for all (x, y)\0 < x < 00, < y < 00 where x — b and y = ^|-. Hence 

A 

- 2 • 



the series representation of I3 is convergent for all a, b, z and <Hor^^±^, A = 0, 1, .... When v is an integer 



or half-integer then the series representation will contain psi functions and logarithmic terms but from the 
structure in (4.4) one can see that the series will be convergent. 

Now let us examine the complicated form coming from I5 or from (2.4). Here we have the factor 



l2°°\v + mi, a, z, p) with mi > 0. For the case p = | wc have from (3.9) and (4.1) 



/' 

Jo 



u+mt -ax -zx 2 , _ x(oo) 



dx = I 2 v + mi, a, z, - 



a -(u+ mi +i) 1 



Writing (4.5) as a G- function, see Mathai 10 , we have 



r(s)r ( s + 1 \r(v + mi + 1 + s) (^-] ds. (4.5) 



V 4 



r (oo) 



u + mi,a,b,-) = 



— (K+mi+l) 

a r<3,o 

■^0,3 







az 2 




4 


0,l,i/+mi+l 



(4.6) 



The behavior of this G-function for small and large values of 2|- is available from Mathai and Saxena or 
Luke 13 . For small values of 2|- this G-function behaves like unity and hence behaves like - ■ 
For large values of ^J- the in (4.6) behaves like 



az 



2\ 6 



az' 



Hence for checking the convergence of I 5 in (2.4) it is sufficient to examine the two series 



oo 2m 

E E 

m— mi— 



2m 

mi 



oc 

E 

m=0 



CO 



2 \ m 



1 - 



ca 



-1 




for 


£(i-i) 




.9 V ca / 



< 1 



9 



or for i — \g\ < c < i + |g| which is equivalent to < c < 1 + \g\ for \g\ > 1, when a = 1, and 



E E 

m— mi— 



C 2 



HI 



mi 



1 

ac 



2\ 3 



= E 



7/ 



m=0 



1 + 7?- 



whcrc 



for 



1 / z 



\2a) 



< 1 or for 



2a[c- \g\]i <z< 2a[c+\gp. 
For example, combined with the condition for small values, we have < z < 2(2 + \g\) for a = 1, |g| > 2. 

For other values of p also the procedure remains the same. We may observe that instead of the I5 in 
(1.9) we can also consider a more general integral of the type 



[( C -y)2+ 5 2]. 



F dy, ^ a, z, c, p, d > 0. 



In this case replace 



[( C -y)2 + ff 2]<*- r(d) 

Then the modification in (2.3) becomes 



1 f 00 

ijr / ^"V^-^+^da;, R(d) > 0. 
(") ./o 



/ x m+d - 1 e- 92a: dx = (. 9 2 )-(" l+d )r(m + d) 
Jo 



(4.7) 



(4.8) 



(4.9) 



and a corresponding expression for (2.4) is obtained. Convergence conditions can be checked exactly the 
same way as in the case of d — 1 = 0. 

Iq can be reduced to a form corresponding to (2.4). Convergence conditions can be checked by converting 
the series form into one of the standard triple series discussed in Srivastava and Karlsson 11 or by using the 
general procedure for a multiple series or by writing the kernel function as a G-function, as described above, 
and then checking the behavior for large and small values of the argument of this G-function. Note that I7 
can be expressed in terms of as follows: 



2 

00 

^ = E(-!) r 

m=0 



2m 



111 



mi— 



/2m 

u 



(-l) mi c 



mi 2m— mi 



^ ~~7 _ 4° C '' { V + m l + ^ fl 5 Z ' P)' 



(4.10) 



n=0 



Writing I^*^ in terms of a G-function and then checking the behavior of the G-function for small and large 
values of the argument one can verify the existence of It. For small values, the G-function in I 2 behaves like 
unity and then for p = \ 

2m 



J 7 = jP(_l)m( rf )^( g 2)-(m+rf) ^ 
m—0 mi— 



2m 
mi 



(-l) mi c 



mi Jlm—vfi\ 



E 



n=0 



>+l) 



e~ a g' 2 \F \d;;~[l- 



ca 
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for ^ [l — i] < 1 or ~ — |<7| < c < i + \g\. For large values the behavior of the G-function is available 
from (5.2) later on and in this case the conditions for the existence are given in (5.8) later on. 

When p = — , m,n = 1, 2, ... it is easy to see that the H-function in all the integrals I\ to ig reduce to 
G-functions, of course with more parameters. General series representations of all forms of G-functions are 
available from Mathai 10 . 



V. Computational aspects 



Anderson, Haubold and Mathai 14 looked into the computational aspects of the integrals I\ to I4. Exact 
computations and graphs are given there for p = i, a = 1 in the cases of I\ for v = 1,2; I2 for v = 0, d = 1, 5; 
I 3 for {u,8,b) — (0,2,0.001), (1,5,1); I4 for {v,t) = (0,1), (0,5). For large values of z one can use 
the asymptotic forms of the integrals for computational purposes. These can be worked out by using the 
asymptotic form of the G-function. From (3.9) and (4.1) 



r(°°) 



v, a, z, — 

' ' ' 2 



-("+1) 



^3,0 

°0,3 



2 

az z 




4 


i/+l,a,i 



(5.1) 



But for large values of 2 4- we have 



G 



3,0 

0.3 



have 




2 




az z 




4 





(5.2) 



By substituting this expression for the G-function in 1\ to Ij we get the following forms for large values of 
^ with a= 1, p=\: 



where rj = 1 — ~ (|) 



for I7I < 1 where 



(!) 



2 / z 

1 



d v+1 e- d 



-Y 

Ad 



-3Mf 



2\ 6 



2\ 3 



-f-V 

o 2 V3/ 



4 



/fi = for d = 1 



v/3 



2\ 2 + 3 



9\ 



-3 



x e 



-MV (d; 



1-ifA 
c V2a 



(5.3) 
(5.4) 
(5.5) 
(5.6) 

(5.7) 



For a broad overview of the aspects of numerical evaluation of special functions, including available software 
packages for this purpose, see on the World Wide Web: http : / /math . nist . gov /nest / . 
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